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1.0  INTRODUCTION 

The  goals  of  this  research  project  were  two-fold.  First,  a  reliable, 
stable  methodology  for  estimating  strong  ground  motion  using  dynamic 
raytracing  in  three  dimensional  geologic  structures  needed  to  be 
developed  and  then  implemented  into  an  interactive  and  efficient 
software  package  at  the  Air  Force  Geophysics  Laboratory.  The  second 
objective  was  to  utilize  that  methodology  to  analyze  the  effects  of 
specific  structural  variations  on  simulated  strong  ground  motion  in 
geologic  basins.  This  report  summarizes  in  detail  the  theoretical  and 
mathematical  basis  for  the  techniques  developed  and  implemented  in  the 
DYNARAY  dynamic  raytracing  software.  In  addition,  this  report 
expands  upon  the  discussion  of  structural  effects  contained  in  a 
previous  report  (Mellman  et  al,  1983)  with  a  discussion  of  the  simulated 
ground  motion  variation  across  a  realistic  Basin  and  Range  structure 
(Yucca  Flats,  NTS). 

The  problem  of  developing  a  methodology  for  estimating  seismic 
amplitudes  in  a  ray  theoretic  approach  which  is  reliable  in  three 
dimensional  varying  structures  was  vastly  more  difficult  than  originally 
believed.  This  effort  required  far  more  research  man-hours  and 
expense  than  estimated  and,  in  the  end,  much  more  than  was  available 
in  the  project  budget.  Completion  required  a  substantial  expenditure  of 
internal  Sierra  R&D  funds  in  addition  to  the  project  budget.  The  end 
product,  however,  is  a  very  sophisticated  solution,  and  potentially  a 
very  useful  software  package.  Much  of  the  difficulties  encountered  by 
the  project  team  resulted  from  the  experience  that,  when  tested 
thoroughly  and  applied  to  realistic  structures,  the  previously  published 
methodologies  for  estimating  ray  amplitudes  failed  at  some  level.  Some 
of  these  problems  were  quite  subtle  and  thus  progress  in  this  area  was 
difficult  and  often  frustrating.  Section  II  of  this  report  discusses  the 
possible  approaches  that  can  be  and  have  been  employed  including  a 
local  WKBJ  methodology  developed  at  Sierra  which,  with  modification, 
eventually  was  adopted  as  the  optimal  technique.  The  DYNARAY 
software  actually  provides  the  user  with  alternative  solutions  allowing 
the  user  to  examine  which  approach  is  most  appropriate  within  the 
constraints  of  model  complexity  and  computational  effort. 
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Section  III  presents  the  results  of  DVNARAY  calculations  for  a  realistic 
basin  model.  The  possible  combinations  and  permutations  of  structural 
variations  are  obviously  enormous  and  in  conjunction  with  the  Project 
Office,  it  was  decided  to  proceed  with  this  part  of  the  project  in  two 
phases.  In  the  first,  a  suite  of  geometrically  simple  basins  were 
defined  and  used  to  help  develop  a  basis  for  understanding  ground 
motions  in  more  complex  structures.  The  modeling  results  for  those 
structures  were  presented  earlier  (Mellman  et  al,  1983).  The  second 
half  of  this  task  was  to  use  a  realistic  model  and  try  to  understand  the 
patterns  of  ground  motion  resulting  from  different  source  locations.  To 
this  end,  we  chose  the  Yucca  Flats  basin  at  the  Nevada  Test  Site  as 
the  ideal  candidate.  This  Basin  is  a  fairly  complex  geologic  structure 
that  has  been  studied  in  detail  by  Herrin  and  his  colleagues  at  SMU. 
The  sediment-Paleozoic  contact  has  been  well  defined  across  the  entire 
basin,  thus  providing  an  excellent  test  case.  As  discussed  in  detail 
later  in  this  report,  the  actual  observed  pattern  of  ground  motion  is 
complex  but  generally  understandable.  Moreover,  this  modeling  effort 
has  identified  a  number  of  important  experimental  considerations  to 
include  in  all  future  modeling  studies  of  this  type. 
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2.0  THEORY 
2.1  MODELS 

To  a  large  degree,  the  type  of  model  permitted  in  a  three-dimensional 
ray  tracing  program  determines  the  efficiency,  accuracy  and  usefulness 
of  such  a  program.  We  wish  the  model  to  be  sufficiently  general  to 
accurately  model  realistic  geologic  structures,  while  minimizing  the 
number  of  grid  points  to  be  stored.  In  addition,  we  wish  to  maintain  a 
representation  of  the  model  that  permits  a  simple  means  of  determining 
the  intersection  points  of  rays  with  layer  boundaries,  since  this  will 
have  a  major  effect  on  program  efficiency.  Finally,  we  wish  the  form  of 
the  model  to  not  introduce  artificial  discontinuities  where  none  exist  in 
the  actual  structure,  since  current  methods  of  calculating  ray 
amplitudes  are  invalid  in  such  models.  Thus  faceted  models,  so  common 
in  finite  difference  techniques,  are  inappropriate  for  ray  methods. 

The  most  general  model  type  would  be  one  composed  of  irregular 
inhomogeneous  layers  with  arbitrary  non-constant  velocity  gradients 
within  each  layer.  For  greatest  generality,  layers  could  be  allowed  to 
pinch  out,  and  thus  not  necessarily  exist  over  the  entire  extent  of  the 
model.  Also,  layers  could  possibly  be  multi-valued,  to  allow  for 
faulting,  overthrusts,  etc. 

The  most  straightforward  realization  of  such  a  model  is  to  specify 
velocity  at  each  point  on  a  three-dimensional  grid.  Such  models  are 
extremely  difficult  to  specify,  since  for  many  realistic  structures 
specification  must  be  virtually  point  by  point  which  is  clearly  not 
practical  in  three-dimensional  models  of  any  substantial  size.  Also,  in 
three  dimensions  the  storage  requirements  for  such  models  are 
enormous,  ruling  out  implementation  on  any  but  virtual  systems  without 
large  overhead  expense.  Further,  ray  tracing  in  such  models  is 
extremely  expensive,  since  it  involves  pointwise  integration  of  a  system 
of  raytracing  equations.  Initial  SGI  efforts  to  produce  a  three 
dimensional  ray  tracing  program  involved  modification  of  a  program 
originally  written  by  Bruce  Julian,  which  does  use  such  models.  It  was 
quickly  evident,  however,  that  this  method  was  much  too  cumbersome 
and  far  too  expensive  to  be  practical  for  modeling  three-dimensional 
basin  structures. 
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In  order  to  obtain  an  efficient,  easily  used  method  the  following 
requirements  were  imposed  on  the  model: 

1)  Homogeneous,  irregular  layers 

2)  Layer  boundaries  are  defined  on  an  even  X,Y  grid  by 
specifying  depth  at  each  grid  point 

3)  Layer  boundaries  are  defined  for  all  X,Y  grid  points 

4)  Layer  boundaries  are  defined  between  grid  points  by  a  9  pt. 
quadric  fit. 

Such  models  offer  a  number  of  advantages.  By  using  layers,  it 
becomes  possible  to  specify  the  model  by  specifying  layer  boundaries 
and  layer  elastic  constants.  This  requires  storage  of  two-dimensional 
rather  than  three-dimensional  arrays,  and  greatly  reduces  storage 
requirements.  Use  of  even  grid  spacing,  with  different  spacing  in  the 
X  and  Y  directions,  and  the  requirement  that  layers  be  defined  for  all 
X  and  Y,  allows  the  use  of  very  efficient  algorithms  to  determine 
ray-boundary  intersections.  The  use  of  smooth  interpolation  between 
grid  points  ensures  reasonable  accuracy  in  amplitude  calculations,  and 
avoids  problems  in  both  ray  capture  and  amplitude  calculations 
introduced  by  artificial  shadow  zones. 

The  model  requirements  above  are  not  as  restrictive  as  they  may  at 
first  appear.  Pinch-outs  may  still  be  realized  by  allowing  two 
boundaries  to  have  the  same  depth  over  a  portion  of  their  range.  It  is 
a  relatively  simple  matter  to  recognize  when  this  occurs,  and  to  compute 
reflection  and  transmission  coefficients  and  ray  paths  appropriate  for 
the  zero  thickness  case.  Similarly,  by  introduction  of  additional  zero 
thickness  layers,  recumbent  structures  may  be  created.  We  also  note 
at  this  time  that  the  constant  velocity  requirement  may  be  relaxed  to  a 
constant,  or  analytic,  velocity  gradient  requirement  without  penalty  in 
storage  but  with  some  penalty  in  execution  time.  This  extra  generality 
was  not  included  in  the  DYNARAY  package. 

The  geologic  model  itself  is  not  a  full  description  of  a  ray  tracing 
model.  We  also  need  to  specify  source  location  and  type,  and  receiver 
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location,  and  ray  type.  DYNARAY  allows  both  explosion  and  double 
couple  sources  at  arbitrary  locations,  and  receivers  at  arbitrary  depths 
and  locations  (although  all  receivers  must  lie  on  a  single  horizontal 
plane).  Both  P  and  S  waves  are  permitted,  with  an  arbitrary  number 
of  reflections  and  conversions  allowed.  Correct  elastic  reflection  and 
transmission  coefficients  are  used  at  all  times.  This  compromise  of 
model  generality  and  program  efficiency  was  felt  optimal  in  providing  a 
practical,  useful  and  accurate  program. 

2.2  KINEMATICS 

In  this  section,  we  consider  the  problem  of  finding  all  ray  paths 
connecting  source  and  receiver  points,  and  the  travel  times  associated 
with  each  ray  path.  In  the  next  section,  we  will  discuss  methods  by 
which  amplitudes  for  each  ray  may  be  determined. 

The  first  problem  that  must  be  considered  is  specification  of  ray  type. 
At  each  interface  encountered  by  a  ray  a  decision  must  be  made 
whether  to  proceed  by  transmission  or  reflection  and  whether  mode 
conversion  will  occur.  To  this  end,  a  ray  description  code  is  used  to 
specify  ray  behavior  at  each  interface.  The  ray  code  used  in 
DYNARAY  has  been  designed  for  notational  compactness.  Transmission 
with  no  mode  conversion  is  assumed  at  each  interface  unless  a  ray 
instruction  explicitly  appears  for  that  interface.  These  instructions  are 
carried  out  sequentially,  as  the  ray  encounters  the  interface  to  which 
they  refer.  Instructions  include  mode  conversion,  transmission  and 
reflection.  Multiple  reflections  from  a  single  layer  may  be  specified. 
The  ray  instruction  code  allows  any  order  multiple  reflection  to  be 
calculated,  with  any  number  of  mode  conversions.  A  more  detailed 
description  of  the  ray  code  may  be  found  in  the  User  Notes. 

Currently,  except  for  direct  arrivals  and  primary  reflections,  ray  code 
generation  must  be  done  manually.  This  is  done  principally  because  of 
the  enormous  number  of  distinct  multiple  reflections  and  conversions 
that  exist  in  three-dimensional  problems.  Unlike  the  flat  layered  case, 
kinematic  and  dynamic  analogs  do  not  exist  for  three-dimensional  ray 
sets  and  automatic  multiple  generation  programs,  which  indiscriminantly 
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generate  all  multiples  of  a  given  order,  tend  to  generate  far  too  many 
ray  codes  to  be  practical . 

Having  generated  a  ray  code,  we  now  consider  the  problem  of 
propagating  a  ray  of  the  appropriate  type  through  the  model.  In  our 
irregular,  homogeneous  layered  model  it  is  a  relatively  simple  task  to 
shoot  a  ray  from  a  specified  source  point  at  a  specified  takeoff  angle 
and  azimuth.  Given  an  initial  ray  direction  t,  we  note  that  in  a 
homogeneous  medium  £  is  constant.  Thus,  t  may  only  be  changed  by 
interaction  with  a  boundary.  At  such  a  boundary,  let  £.  be  the  ray 
vector  for  the  incident  wave  at  the  boundary,  f  the  outgoing  ray 
vector,  and  h  the  surface  normal  at  the  intersection  point  with 
fj  •  n  <  0  and  c.  and  cq  the  incident  and  outgoing  material  velocities. 
Then 

t  =  co  fix  (t.  x  fi)  ±  (1  -  co  |  t.  x  fi|^)^  fi  (1) 

c .  2 

1  c . 

i 

if  |  t.  x  fi  |  <  1  and 

c . 

x 

i  =  1  fix  (t.  x  fi)  (2) 

O  |~T - *T  1. 

I  ti  x  n| 

if  !l°  |  ti  X  fi  |  >  1 

c. 

1 

defines  the  outgoing  ray  vector  for  transmission  (reflection).  This  may 
be  recognized  as  a  vector  statement  of  Snell's  law. 

The  problem  of  ray  propagation  may  thus  be  seen  to  reduce  to 
determining  the  successive  intersections  of  lines  with  irregular 
bjundaries,  as  determined  by  the  ray  code.  A  search  algorithm  is  first 
use<-<  o  restrict  the  intersection  range  to  a  single  grid  spacing.  Then, 
siru_e  the  boundary  is  locally  described  by  a  quadric,  an  exact  inter- 
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section  point  may  be  obtained  analytically.  Determination  of  the  surface 
normal  at  this  point  is  also  analytic. 

It  is  possible  to  make  some  simple  extensions  to  the  theory  for  the 
homogeneous  model  although,  as  noted  below,  at  some  cost.  In 
particular,  Hubral  (1978)  discusses  layered  media  with  constant  velocity 
gradients.  Ray  paths  in  such  media  are  circular  arcs  and  the  ray 
tracing  problem  reduces  to  one  of  determining  the  intersection  of  circles 
with  irregular  surfaces.  This  results  in  a  more  complicated  and  time 
consuming  intersection  algorithm.  More  complicated  models,  such  as 
those  discussed  by  Hron  and  Cerveny  (1980),  generally  require  finding 
the  intersections  of  circular  arcs  with  a  number  of  model  cells  of 
constant  velocity  gradient.  Such  methods  are  prohibitively  expensive 
for  use  in  large  scale,  realistic  three  dimensional  problems. 

Having  solved  the  problem  of  propagating  rays  through  the  model,  we 
now  wish  to  address  the  two-point,  or  ray  capture,  problem.  That  is, 
we  wish  to  find  all  rays  for  a  given  ray  code  which  connect  the  source 
and  receiver.  Due  to  the  three  dimensional  nature  of  the  model,  we 
cannot  in  general  even  vaguely  predict  the  take-off  angle,  azimuth  or 
number  of  such  rays.  It  is  therefore  necessary,  no  matter  what 
capture  algorithm  is  used,  to  sample  ray  parameter  space  to  at  least 
obtain  starting  models  for  the  capture  algorithm.  We  do  not,  however, 
want  to  resort  to  saturation  shooting  in  three  dimensional  problems. 

At  this  point,  it  is  of  some  interest  to  note  that  we  very  seldom  wish  to 
solve  the  problem  usually  posed  in  mathematical  treatments  of  ray 
capture.  That  is,  we  seldom  have  problems  involving  a  single  source 
and  a  single  receiver.  Instead,  we  usually  have  problems  involving  a 
single  source  and  many  receivers.  Indeed,  DYNARAY  is  intended  to 
generate  simultaneously  the  response  on  a  large,  regularly  spaced  grid 
in  order  to  produce  contours  of  maximum  ground  response.  Algorithms 
that  are  efficient  for  one  problem  are  often  not  efficient  for  the  other. 
For  the  multiple  receiver  case,  we  note  that  the  same  set  of  rays  used 
to  sample  ray  parameter  space  may  be  used  for  all  receivers.  We  will 
call  this  the  working  ray  set  and  use  this  set  of  rays  to  construct  an 
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efficient  capture  algorithm.  We  note  that  as  the  number  of  receivers 
increases,  the  overhead  associated  with  generating  the  working  ray  set 
proportionately  decreases,  and  that  in  very  large  problems  fairly  dense 
working  ray  sets  can  be  used  with  little  additional  cost  in  efficiency. 

We  may  identify  ray  captures  using  the  working  ray  set  in  the  following 
manner.  We  identify  rays  in  the  working  ray  set  by  their  take-off  angle 
and  azimuth,  (<j>,  0)  and  assume  without  loss  of  generality  that  working 
rays  are  generated  by  incrementing  <{>  and  0.  Let  (q^,  0^)  be  a 

working  ray  which  emerges  at  point  .  If  6q>  and  60  are  increments  of 
angle,  let  (q>2,  02)  =  (q>1  +  6^,  0^  and  (q>3>  03>  =  ( ,  01  +  60)  with 
emergence  points  X 2  and  X^.  We  define  a  capture  if  the  receiver  point 
Xr  is  contained  in  the  triangle  formed  by  X^  X2  and  X3.  Thus,  the 
problem  of  identifying  captures  is  reduced  to  a  check  of  triangles 
defined  by  the  emergence  points  of  the  working  rays.  The  capture 
process  may  then  continue  by  refining  the  ray  parameter  estimate  using 
either  ray  shooting  or  ray  bending  methods.  This  refinement  process 
is  complicated  somewhat  by  the  non-linearity  of  the  problem.  The  true 
captured  ray  may  lie  outside  the  triangle  in  ray  parameter  space  that 
defines  the  capture.  There  may  be  multiple  captured  rays  within  the 
triangle,  or  a  shadow  zone  may  exist  so  that  there  is  no  true  captured 
ray,  although  in  this  case  a  diffracted  arrival  will  still  exist.  The 
capture  process  is  considered  complete  when  a  ray  (<J>p  0^)  is  found 
such  that 


I  Xf  -  Xr  |  <  e. 

We  have  found  that  a  combination  of  search  and  gradient  methods  used 
with  a  ray  shooting  technique  works  reasonably  well  as  a  refinement 
technique.  Any  refinement  technique,  however,  will  be  expensive,  and 
we  would  like  to  avoid  using  such  techniques  if  possible. 

It  is  possible  to  avoid  refining  ray  capture  in  cases  where  the  principle 
interest  is  in  arrival  time  and  amplitude,  rather  than  exact  ray  path, 
and  if  an  amplitude  method  is  used  that  is  not  sensitive  to  very  small 
scale  model  features.  We  will  discuss  such  amplitude  methods  in  the 


SGI-R-83-096 


9 


next  section.  We  will  call  our  alternative  to  refining  ray  capture  "time 
capture".  We  note  that  for  a  ray  (0q,  0  ),  we  can  predict  the  arrival 
time  of  the  true  captured  ray  ((jy,  0^)  at  Xr  by 


(3) 


where  tQ  is  the  travel  time  of  (0o,  0  )  and  pQ  is  the  projection  of  the 
slowness  vector  at  XQ  onto  the  surface.  If  the  rays  defining  the 
capture  all  predict  arrival  times  at  Xr  within  some  error  limit,  we  may 
define  a  time  capture  to  have  occurred.  The  travel  time  is  then  taken 
as  a  weighted  average  of  predicted  times.  Otherwise,  a  refinement 
procedure  may  be  invoked  until  an  acceptable  time  capture  is  obtained. 
In  many  applications,  this  refinement  is  unnecessary.  In  the  present 
application,  where  the  primary  aim  is  to  predict  structural  amplification, 
small  time  errors  are  unimportant.  DYNARAY  uses  the  time  capture 
method,  and  achieves  a  substantial  savings  in  computer  time  compared 
to  more  standard  capture  methods. 

2.3  AMPLITUDES 

In  this  section  we  consider  the  problem  of  determining  the  amplitudes  of 
rays  for  a  given  instruction  set.  Each  amplitude,  together  with  the 
source  time  function  and  information  on  the  direction  of  motion  for  that 
ray,  forms  a  single  ray  arrival.  The  seismogram  is  then  the  sum  of  all 
such  arrivals. 


It  is  not  our  intention  to  provide  a  complete  theoretical  development  of 
all  methods  used  to  determine  amplitudes,  since  such  developments  are 
readily  available  for  most  of  these  methods  in  the  literature.  Rather, 
we  will  concentrate  on  the  physical  basis  of  these  methods,  and  rely  on 
final  expressions  for  ray  amplitudes  taken  from  the  literature  and 
presented  without  formal  justification.  Fortunately,  research  in  ray 
amplitudes  in  three  dimensions  has  been  quite  active  lately,  and 
significant  progress  on  some  classic  problems  in  ray  theory  has 
occurred  and  has  been  presented  in  the  geophysical  literature. 
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Traditionally,  ray  techniques  have  used  amplitudes  derived  from 
geometric  optics.  Such  amplitudes  are  valid  in  illuminated  regions 
where  changes  in  wavefront  curvature  are  small  compared  to  the 
curvature  itself.  Derivations  of  the  equations  governing  ray  tracing 
and  geometric  optics  are  readily  available  in  texts  such  as  Aki  and 
Richards  (1981).  A  more  comprehensive  treatment  is  given  by  Hron 
(1982). 


The  geometric  optics  amplitude  for  a  ray  connecting  source  and  receiver 
is  given  by  Chapman  and  Drummand  (1982)  as 


A(£o,t)  =  Im 


M(£0,t) 

5/2  ±  ±li 

4  na  a  p  p 


R 


6(t-T(Eo)) 


(4) 


where  M(d  ,t)  =  m(po,t)+iH(m(go,t))  is  an  analytic  source  term,  with 
mo(pQ,t)  containing  radiation  pattern  and  source  time  function  terms,  H 
is  a  Hilbert  transform,  the  subscript  o  refers  to  source  values,  a  is 
velocity,  p  is  density,  p  is  the  horizontal  projection  of  the  slowness 
vector,  T(pq)  is  the  travel  time  for  a  ray  with  inital  ray  parameter  pQ 
and  R  is  a  transmission-reflection  product.  The  parameter  o  is  the 
KMAH  parameter  defined  by  Chapman  and  Drummond  (1982),  and 
measures  the  number  of  caustics  encountered  by  the  ray. 


Aside  from  constant  terms,  which  depend  on  the  elastic  properties  of 
the  source  and  receiver  regions,  there  are  three  terms  which  determine 
the  amplitude  of  a  geometric  ray.  These  are: 

1)  Source  radiation  pattern  and  time  function 

2)  Transmission-reflection  product 

3)  Geometric  spreading 

The  first  two  terms  are  well-behaved,  and  are  common  to  all 
ray-theoretic  methods  that  we  will  deal  with  in  this  section.  It  is  the 
geometric  spreading  term  which  causes  the  optics  solution  to  become 
inaccurate  at  caustics  and  is,  in  general,  responsible  for  the  failure  of 
optics  solutions  in  any  situation  where  diffraction  effects  become  impor- 


L 
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tant,  We  will  deal  with  several  methods  of  generalizing  geometric 
spreading  factors  to  avoid  some  of  these  problems. 

First,  however,  we  need  to  address  the  source  and  reflection  terms. 
In  DYNARAY,  two  source  types  are  permitted,  earthquake  (double 
couple)  sources  and  explosion  sources.  The  amplitude,  due  to  radiation 
pattern  effects,  of  the  source  term  may  be  obtained  using  a  saddle 
point,  or  first  motion  approximation  on  the  whole  space  response  to  a 
point  explosion  or  shear  dislocation.  Following  Langston  and 
Helmberger  (1975),  we  find  the  amplitude  contribution  of  the  source 
radiation  pattern  to  be 

3 

A  =  I  6. A.  (0,A,6)  C.  (5) 

j=0  J  3  J 

with  0,A  and  6  the  strike,  rake  and  dip  of  the  fault,  j  =  0  the 
explosion  source,  and 

1  for  j  =  0  and  an  explosion  source 
6=0  for  j  =  0  and  a  double  couple  source 
J  1  for  j  >  0  and  a  double  couple  source 
0  for  j  >  0  and  an  explosion  source 

and 

C0  =  '1/0o 

C1  =  "Po 

C2  =  2epna 
2  2 
C3  =  Po'2na 

Aq  =  1 

A^  =  sin  26  cos  A  sin  6  +  %  cos  20  sin  A  sin  25 
A2  =  cos  0  cos  A  cos  6  -  sin  9  sin  A  cos  26 
A^  =  *4  sin  A  sin  26 
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with 


1  for  downgoing  waves 

£  = 

-1  for  upgoing  waves 


and 

na  =  C1/®2  -  p2)*4 

For  initial  SV  wave  amplitudes  due  to  the  source  radiation  pattern,  we 
have 


A  =  I  6. A.  (0 ,A,6)  SV 

j=l  J  J  J 


(6) 


with 


SV1  =  eprip 
sv2  =  n2  -  P2 


SV3  =  3ep  Hp 
and  np  =  (VP*  -  P2)* 


For  initial  SH  waves,  we  have 


A  =  I  6.  SH  A  3  (0,A,6) 
j=l  J  J  J 


(7) 


rith 


S“l  =  f? 


SH,  =  :i!l 

^  a 


A,  =  cos  0  cos  A  sin  6  ■  *4  sin  20  sin  A  sin  26 
A 

A^  =  -sin  0  cos  A  cos  6  -  cos  0  sin  A  cos  26 
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In  addition  to  the  source  amplitude  term,  the  source  term  m(t,p) 
includes  a  far-field  time  function.  This  time  function  is  the  same  for 
all  rays,  and  is  included  in  the  seismograms  in  a  final  convolutional 
step.  For  simplicity,  we  introduce  the  analytic  time  function 
M(t,p)  =  m(t,p)  +  iH(m),  in  order  to  facilitate  implementation  of  n/2 
phase  shifts  contained  in  both  reflection  and  geometric  spreading  terms. 
By  including  these  as  imaginary  portions  of  the  amplitude,  a  final 
seismogram  containing  correct  phase  shifts  may  be  efficiently  obtained. 

The  reflection  coefficient  term  R  in  equation  4  represents  a  product  of 
all  transmission  and  reflection  coefficients  for  interfaces  encountered  by 
the  ray.  It  may  be  shown  (Popov  and  Psencik,  1978)  that  the  usual 
generalized  plane  wave  reflection  and  transmission  coefficients 
appropriate  for  flat  layered  problems  are  also  appropriate  for  the  first 
order  asymptotic  solutions  in  laterally-varying  media.  The  form  used  in 
DYNARAY  for  these  generalized  reflection  and  transmission  coefficients 
is  that  of  Helmberger  (1968),  and  is  appropriate  for  both  fluid  and 
solid  elastic  media.  As  the  expressions  for  these  coefficients  are 
somewhat  lengthy  and  are  readily  available  in  the  literature,  they  will 
not  be  reproduced  here. 

The  main  problem  encountered  in  applying  the  plane  wave  reflection  and 
transmission  coefficients  is  the  fact  that  neither  ray  parameter  nor  local 
SH  and  SV  directions  are  conserved  in  laterally-varying  media.  In 
irregular,  homogeneously  layered  media,  changes  in  ray  parameter  and 
S  wave  polarization  occur  only  at  layer  boundaries.  It  is  therefore  only 
necessary  to  change  coordinate  systems  at  interfaces  in  this  type  of 
medium  to  be  able  to  correctly  handle  mode  conversions,  reflections  and 
transmissions.  The  problem  is  thus  transformed  into  one  of  successive 
interactions  of  a  ray  with  planes  tangent  to  interfaces  at  the 
intersection  points,  with  local  ray  parameter  and  shear  wave  polarization 
directions  defined. 

It  is  in  general  simplest  to  formulate  reflection  and  transmission 
problems  in  terms  of  ray  centered  coordinates  (Hubral,  1979).  To 
discuss  this  approach,  we  need  to  define  a  few  terms.  First  let  e^  be 


SGI -R -83-096 


14 


the  direction  of  propogation  of  a  ray,  and  and  §2  be  orthogonal 
directions  which  form  the  "natural"  shear  wave  polarization  directions. 
Initially,  we  choose  e^  horizontal  and  specify  that  the  coordinate  system 
is  right-handed.  At  any  interface,  we  perform  a  coordinate  rotation  to 
new  coordinates  (ej,  ej,  e^,)  such  that 


~  f  -  £  % » 

1 

-j  e3  x  fi 

1  —  ft  1  u 

(8) 

3  3  ’  1 

Je3  x  n 

e2  ‘  e3  x 

el 

The  directions 

ej  and 

ei  are 

the  local  SH 

and 

SV  directions 

respectively,  and 
p1  =  ^(1-(h  •  e3 

the  local 
,)^)^  with 

ray  parameter  is  given 
c  the  appropriate  P 

by 

or  S 

wave  velocity. 

Incident  S  waves 

are  first 

resolved 

into  local  SH  and  SV 

by 

ASV  =(e2 

Asv  +< 

>,  • 

ash 

(9) 

a'sh  =(ei 

e'l)  ftSH 

:v*i> 

ASV 

(10) 

where  A^v  is  the  previous  local  SV  amplitude  (which  may  be  complex) 
and  Ash  is  the  previous  local  SH  amplitude.  The  appropriate  reflection 
or  transmission  coefficients  are  then  applied,  to  obtain  new  P  or  S 
amplitudes,  with  outgoing  direction  e^  given  by  equation  1  and 

e”  =  e[  e'2  =  x  e”  (11) 

This  new  cooordinate  system  is  maintained  until  the  next  layer  is 
encountered,  at  which  time  the  procedure  is  repeated.  We  note  that  in 
media  with  velocity  gradients,  the  e^  and  e ^  axes  rotate  about  e^  as  the 
wave  propagates.  This  is  one  reason  that  ray  tracing  in  a  general 
inhomogeneous  medium  is  more  expensive  than  in  our  restricted  medium, 
since  this  rotation  must  in  general  be  found  by  integration  along  the 
entire  ray  path. 
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For  ray  sets  involving  shear  phases,  it  is  necessary  to  either  maintain 
complex  SV  and  SH  amplitudes  or  to  specify  a  single  complex  S 
amplitude  and  appropriate  polarization  information.  Shear  wave 
polarizations  will  be  linear  unless  a  critical  angle  is  encountered,  at 
which  time  elliptic  polarization  is  possible.  In  any  case,  travel  time  and 
geometric  spreading  will  be  identical  for  SV  and  SH  as  long  as  the 
medium  is  isotropic.  A  discussion  of  anisotropic  media  is  beyond  the 
scope  of  this  report. 

In  general,  the  transmission-reflection  product  will  be  complex. 
Transmission  coefficients  will  themselves  be  real  unless  a  critical  angle 
is  encountered.  In  our  computational  approach,  when  a  critical  angle  is 
encountered  on  transmission,  the  ray  is  simply  terminated  and  zero 
amplitude  assigned.  Post-critical  reflections  are  permitted,  and  in  this 
case  the  reflection  coefficient  will  become  complex.  Physically,  this 
complex  coefficient  corresponds  to  a  phase  shift  and  may  be  represented 
by  using  a  weighted  sum  of  the  source  wavelet  and  its  Hilbert 
transform.  This  is  accomplished  using  the  analytic  source  term  in 
equation  4,  as  previously  described.  The  result  is  an  asymptotically 
correct  result  for  post-critical  reflections,  similar  to  that  obtained  for 
flat-layered  media  using  a  first  motion  approximation  about  reflection 
time  in  a  Cagniard-type  solution. 

The  third  factor  in  equation  4  influencing  amplitudes  is  the  geometric 
spreading  factor,  L.  In  the  optical  solution,  we  have 


‘*5 


L  = 


dx 

I  t  •  i  | 

d*o 

(12) 


The  phase  factor  o(x,xo),  called  by  Chapman  the  KMAH  index,  is  the 
number  of  caustics  encountered  by  the  ray.  Each  caustic  introduces  an 
additional  n/ 2  phase  shift,  with  a  point  caustic  introducing  a  n  phase 
shift.  The  identification  of  caustics  is  rather  difficult  in  three 
dimensions,  and  is  best  handled  by  examining  changes  in  the  wavefront 
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curvature  matrix  along  the  ray.  We  will  discuss  that  approach  in  our 
discussion  of  dynamic  ray  tracing  later  in  this  report. 

The  determinant  j  d  x/d  g^  |  in  the  spreading  factor  is  ordinarily 
determined  by  shooting  additional  rays  at  small  increments  of  ray 
parameter  and  differencing  emergence  points  to  estimate  the  necessary 
derivatives.  This  is  the  method  used  by  Hong  and  Helmberger  (1978) 
in  the  method  they  called  Glorified  Optics.  An  alternative  to  this 
method  of  evaluating  the  determinant  was  proposed  by  Hubral  (1978). 
This  method,  usually  called  dynamic  ray  tracing,  makes  use  of  the 
wavefront  curvature  matrix,  K,  expressed  in  ray  centered  coordinates, 
to  determine  the  geometric  spreading.  Yet  another  alternative  and  the 
one  we  will  use  is  due  to  Cerveny  (1983). 

The  matrix  K  is  a  2  x  2  matrix  which  gives  the  relative  location  along 
the  e^  and  e~  directions  of  rays  at  nearby  ray  parameters 
g  +  and  g  +  6g,,.  The  inverse  of  K,  which  we  will  call  R,  is  the 
radius  of  curvature  matrix.  Cerveny  (1983)  demonstrated  that  if  we  let 
K  =  PQ  1 ,  then  P  and  Q  satisfy 


39  = 

dS 


cP 


(13) 


where  V  is  a  matrix  describing  the  second  derivatives  of  the  velocity 
function  perpendicular  to  the  propagation  direction,  S  is  the  arc  length 
along  the  ray  path  and  all  matrices  are  in  ray-centered  coordinates. 
For  irregular  homogeneous  medium,  V  =  0  and 


dP 

dS 


0 


39  = 

dS 


cP 


(14) 
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Physically,  the  inverse  of  the  P  matrix  describes  the  change  in  ray 
parameter  in  a  ray  tube  about  the  central  ray  for  a  change  in  initial 
ray  parameter.  The  matrix  Q  is  the  change  in  spatial  location  of  a  ray 
for  a  change  in  initial  ray  parameter.  The  determinant  of  Q  may  be 
related  to  the  determinant  of  dx/dg^  by 


dx 

I  £  *  2  I  =  det  Q 
Thus  the  spreading  factor  L  is  given  by 


(15) 


(det  Q)'i 


(16) 


The  initial  values  of  P  and  Q  are 

P  =  I 
Q  =  O 

where  I  and  O  are  the  2x2  identity  and  zero  matrix 
then  be  determined  by  integration  of  equation  14. 

In  a  homogeneous  layer  this  integration  is  quite  simple  as  P  is  constant 
and  Q  is  thus  given  by  the  value  of  Q  on  entering  the  layer  plus  some 
constant  multiple  of  P.  At  boundaries,  the  transformation  of  P  and  Q 
may  be  obtained  by  matching  phase  and  displacement  across  the 
boundary.  The  expressions  for  these  transformations  are  given  by 
Hubral  (1980)  and,  as  they  are  somewhat  lengthy,  will  not  be  repeated 
here. 

The  use  of  the  matrix  Q  also  provides  a  convenient  method  of 
evaluating  the  KMAH  factor  o(x,xQ),  since  caustics  correspond  to  zero 
eigenvalues  of  Q.  In  a  homogeneous  medium,  it  is  only  necessary  to 
examine  the  changes  in  sign  of  the  determinant  and  trace  of  Q  as  the 
ray  enters  and  leaves  each  layer  to  get  a  correct  value  of  a.  This 
provides  a  very  efficient  means  of  getting  the  correct  phase  shift  even 
in  very  complicated  structures. 


(17) 

P  and  Q  may 
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Dynamic  ray  tracing,  as  described  here,  is  one  of  the  options  for 
determining  amplitudes  provided  in  DYNARAY.  However,  as  the 
solution  is  not  valid  at  or  near  caustics  or  in  shadow  zones,  several 
other  methods  are  also  provided.  One  such  method,  which  is  valid  at 
caustics,  is  the  WKBJ  seismogram  initially  described  by  Fraser  and 
Phinney  (1980)  and  Sinton  and  Fraser  (1982)  and  generalized  to  a  fully 
uniform  theory  using  Maslov  asymptotic  theory  by  Chapman  and 
Drummond  (1982). 

Again,  it  is  not  our  intention  to  give  a  full  theoretic  treatment  of  WKBJ 
seismograms,  as  this  has  already  been  done  at  great  length  by  Chapman 
and  Drummond  (1982).  Instead,  we  will  briefly  discuss  the  physical 
meaning  of  the  WKBJ  solution  and  give  the  final  expression,  drawn  from 
the  literature,  used  in  DYNARAY. 


The  WKBJ  seismogram  uses  a  decomposition  of  a  point  source  into  plane 
wave  components.  Each  plane  wave  component,  for  a  given  ray 
instruction,  is  associated  with  a  ray  of  appropriate  take-off  angle.  The 
final  displacement  is  built  up  as  an  integral  over  initial  ray  parameter, 
with  each  ray  parameter  contributing  equally  at  a  time  determined  by 
the  intersection  of  a  plane  perpendicular  to  the  ray  propagation 
direction  with  the  receiver  point. 


Mathematically,  the  WKBJ  response  for  a  P-wave  is  given  by 


u  =  Ilm 
ray 
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where  t  is  a  reduced  travel  time  for  each  ray  with 


t  =  T-  p . (x-x  ) 
*  -  -o 


The  expression  for  shear  waves  is  similar,  but  with  appropriate  changes 
in  the  direction  of  motion  and  elastic  parameters.  A  discrete  realization 
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in  ray  parameter  space  of  equation  18  is  easily  evaluated.  This  gives 
us 


Re  2  °  I!  5£0  f  B(t) 

u  S  I  Ira  I  —  *  -  (191 

1%,  ®» dt  »"2  p V C'2  i  £  •  *  i*  1 1-  p  2 il 

where  B(t)  is  a  boxcar  function  of  unit  area  extending  in  time  from  the 
minimum  to  the  maximum  reduced  time  in  each  discrete  dp^  cell. 

We  note  that  when  using  WKBJ  as  given  in  equation  19,  ray  capture  is 
no  longer  necessary.  Rather,  the  response  is  built  up  by  summation 
over  the  working  ray  set.  This  is  somewhat  deceptive,  however,  for 
unless  fine  binning  is  used  in  ray  parameter  space,  travel  time  errors 
can  result.  This  error  is  comparable  to  the  error  in  the  first  order 
time  capture  method.  Further,  fine  binning  is  required  in  order  to 
obtain  a  smooth  result  from  equation  19.  This  is  the  penalty  paid  to 
retain  the  full  frequency  dependence  of  the  WKBJ  solution.  In  general, 
the  fine  binning  requirement  makes  the  WKBJ  seismogram  considerably 
more  expensive  than  the  dynamic  ray  tracing  seismogram,  even  for 
complicated  models  with  multipathed  arrivals. 

The  range  of  validity  of  the  WKBJ  seismogram  has  been  extensively 
discussed  by  Sinton  and  Fraser  (1982)  and  Chapman  and  Drummond 
(1982).  WKBJ  seismograms  are  valid  at  caustics  and  agree  with  optics 
solutions  when  both  are  valid.  WKBJ  becomes  inaccurate  when  dp/dp^ 
is  zero  or  infinite.  This  occurs  at  so-called  "telescopic"  points,  where 
rays  are  parallel,  and  corresponds  to  a  caustic  in  ray  parameter  space. 
In  addition,  WKBJ  produces  plausible,  although  not  necessarily  correct, 
results  in  shadow  zones.  Shadow  zone  arrivals  will  not  generally  have 
the  correct  amplitude  and  may  have  incorrect  travel  times  as  well,  since 
the  travel  time  information  is  based  on  extrapolation  away  from  rays 
through  unsampled  regions  of  the  model.  We  note  that  not  only  shadow 
zones,  but  any  discontinuity  in  reduced  travel  time  will  produce  such 
arrivals.  These  will  generally  be  small  and  of  little  consequence  for  the 
strong  ground  motion  assessment  problem  of  interest  to  us  here,  but 
can  cause  serious  problems  in  exploration  applications. 
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Despite  the  cost  and  the  spurious  arrivals,  WKBJ's  behavior  in  caustics 
makes  it  a  valuable  technique,  and  it  was  included  in  DVNARAY  as  an 
option.  It  is  reasonable  to  ask,  however,  whether  some  less  expensive 
local  approximation  to  WKBJ  exists  which  is  comparable  to  optics  or 
dynamic  ray  tracing  in  its  use  but  does  net  fail  catastrophically  at 
caustics.  We  will  call  such  an  approximation  a  "local  WKBJ 
approximation" . 

The  key  to  the  local  WKBJ  approximation  lies  in  the  fact  that,  for  a 

minimum  or  maximum  time  geometric  arrival  away  from  a  caustic,  the 

amplitude  of  the  arrival  is  determined  by  the  ray  parameters  in  the 
neighborhood  of  the  geometric  ray  parameter.  Other  ray  parameters  in 
the  integral  merely  serve  to  prevent  the  amplitude  from  changing.  Let 

6  g£  define  a  closed  region  about  g^  such  that  x  =  t  t  e  on  the 

boundary  of  the  region  for  a  minimum  (maximum)  time  phase.  Then  the 
ray  amplitude  associated  with  that  phase  is  given  by 


■in 


A(w)  = 


(meas  (6p  )) 


2  1,  \  L  5  / 

8 n  p1  p  1  a1  a  ' 
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2  (t  •  z)*(l-  p  2)* 


The  region  Sg^  may  be  estimated  by  fitting  a  quadratic  curve  in  6p^ , 
6p^  to  the  x  surface.  This  may  be  done  either  using  adjacent  working 
rays,  or  a  new  group  of  rays  shot  at  some  pre-specified  increments  of 
take-off  angle.  If  the  increments  of  take-off  angle  and  azimuth  are 
chosen  such  that  the  resulting  emergence  points  lie  along  the  principle 
curvature  directions,  a  particularly  simple  form  results  for  the 
quadratic  fit,  in  that  the  cross  terms  will  be  zero.  The  region  of  6g£ 
will  then  be  an  ellipse  with  axis  lengths  given  by  the  quadratic 
coefficients  times  e,  and  an  area  of  nz2  times  the  product  of  the 
coefficients . 

In  cases  where  the  x  surface  contains  a  saddle,  we  determine  the 
amplitude  of  the  Hilbert  transform,  rather  than  the  amplitude  of  the 
true  arrival.  In  this  case,  the  coefficients  of  our  fit  to  the  x  surface 


| 
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will  have  opposite  signs.  By  changing  the  sign  of  the  positive 
coefficient  and  incrementing  o  by  one,  we  in  effect  determine  the 
amplitude  of  the  Hilbert  transform  of  the  true  arrival,  without  adding  to 
the  complication  of  our  algorithm. 

In  the  neighborhood  of  a  caustic,  the  t  surface  will  no  longer  be  locally 
quadratic ,  and  the  amplitude  in  equation  20  will  depend  on  l.  Equation 
20  is  still  a  valid  approximation,  but  the  amplitude  A  will  be  an  average 
amplitude  over  a  time  interval  of  duration  e.  We  are  thus  replacing  a 
frequency  dependent  amplitude  with  an  average  constant  amplitude, 
which  will  be  quite  adequate  for  most  purposes.  This  will  be 
particularly  true  if  the  source  wavelet  is  narrow  band  and  e  is  chosen 
as  the  inverse  of  the  predominant  period. 

As  a  quadratic  fit  to  the  t  surface  is  not  valid  in  the  neighborhood  of  a 
caustic,  some  alternative  method  must  be  found  for  estimating  meas  (6 
jd£).  The  most  straightforward  way  is  to  fit  a  cubic  to  the  t  surface. 
This,  however,  results  in  a  rather  awkward  expression  for  An 

alternative  is  to  do  piecewise  quadratic  fits.  At  a  caustic  point,  four 
such  fits  will  be  necessary,  with  each  valid  in  a  sector.  The  character 
of  the  t  surface  will  usually  change  from  minimum  or  maximum  time  to 
minimax  behavior  between  these  sectors.  The  amplitude  contribution 
and  phase  shifts  of  each  sector  must  be  evaluated  separately.  It  is  of 
some  interest  to  note  that  if  the  piecewise  quadratic  approximation  is 
used,  it  is  possible  to  rewrite  equation  20  as 


-in 

1  ° 

A(U))  =  - * - - -  I  — - ^  (21) 

4n  p*  pj*  ctj2  (f.  •  £)*  sectors 

where  6X£  is  the  region  swept  out  by  emergent  rays  in  a  time  ±e  of  tQ. 
This  expression  is  a  very  close  analog  of  equation  4  which  is 
well-behaved  in  the  neighborhood  of  caustics. 

The  local  WKBJ  amplitudes  are  offered  as  the  default  amplitude  method, 
together  with  time  capture,  in  the  DVNARAY  solution.  While  losing  the 
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inherent  frequency  dependence  of  the  full  WKBJ,  the  local 
approximation  is  much  faster,  reasonably  accurate,  and  avoids  problems 
of  spurious  arrivals.  While  not  implemented  in  the  current  version,  it 
would  be  possible  to  obtain  an  approximation  for  diffracted  arrivals  from 
local  WKBJ  by  expanding  all  extremal  points  on  the  i  surface,  not  just 
those  corresponding  to  geometric  arrivals.  This,  however,  remains  a 
topic  for  future  research. 

2.4  RESPONSE  SPECTRA 

The  goal  of  the  mathematically  complex  operations  arising  from  the 
theory  summarized  in  the  preceding  portions  of  this  report  is  to  provide 
the  basic  input  to  estimates  of  lateral  variations  in  strong  ground 
motion.  DVNARAY  utilizes  the  methodologies  summarized  here  to 
produce  synthetic  accelerograms  for  either  earthquake  or  explosive 
sources.  However,  these  individual  accelerograms  need  to  be  combined 
and  reduced  to  a  format  useful  in  engineering  estimates  of  strong 
ground  motion  hazards.  The  format  normally  used  for  that  purpose  is 
the  pseudo  velocity  response  spectra.  The  uses  and  computation  of 
response  spectra  are  well  described  in  numerous  published  books  and 
articles  on  earthquake  engineering.  For  the  benefit  of  the  reader,  a 
very  brief  description  of  response  spectra  is  presented  below. 

Pseudo  velocity  response  spectra  are  developed  through  the  computation 
of  the  maximum  displacement,  relative  to  the  support  structure,  of  a 
damped  simple  harmonic  oscillator  subject  to  the  specified  strong  ground 
motion  time  history.  Each  period  represented  in  the  response  spectrum 
is  derived  from  a  separate  analysis  of  an  oscillator  with  the  same 
natural  period.  Because  many  structures  can  be  crudely  approximated 
by  a  system  of  damped  simple  harmonic  oscillators,  response  spectra  are 
a  simple  but  useful  analysis  of  strong  ground  motion  data.  At  high 
frequencies,  or  equivalently  small  periods,  it  is  useful  to  note  that  the 
motions  of  the  oscillator  tends  towards  the  values  of  the  input  ground 
motion.  In  the  limit  of  zero  period  the  oscillator  becomes  perfectly 
rigid  and  the  motions  of  the  oscillator  exactly  reproduce  the  input 
ground  motion.  Hence  the  terms  peak  acceleration  and  zero  period 
response  are  often  used  interchangeably .  For  sinusoidal  motion  the 
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difference  between  the  velocity  and  the  displacement  of  an  oscillator  is 
simply  a  factor  of  the  angular  frequency  ui.  The  pseudo  velocity 
response  is  computed  by  multiplying  the  peak  relative  displacement  of 
the  oscillator  by  its  natural  frequency.  It  is  well  recognized  that  for 
complex  motions  this  approximation  can  deviate  from  the  true  maximum 
velocity  response  of  the  oscillator. 
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3.0  APPLICATION  OF  DYNARAY  TO  A  GEOLOGIC  BASIN 
3.1  INTRODUCTION 

Estimation  of  strong  ground  motion  parameters  for  engineering  design 
(eg:  time  histories,  response  spectra,  strain  and  so  on)  generally 
require  a  best  mean  estimate  and  a  measure  of  the  dispersion  or 
uncertainty  of  the  mean.  For  critical  facilities,  such  as  nuclear  power 
plants,  design  is  often  evaluated  on  the  basis  of  expected  motions  that 
exceed  the  mean  by  one  standard  deviation  (i.e.,  84%  of  the  data  fall 
below  this  level).  For  source  locations  and  magnitudes  that  can  be 
estimated  in  advance,  such  as  the  maximum  magnitude  earthquake 
located  along  the  closest  approach  of  a  fault,  and  for  site  conditions 
that  are  carefully  quantified,  common  estimates  of  dispersion  range  from 
factors  of  1.5  to  2.0.  This  multiplicative  factor  essentially  represents 
inherent  uncertainties  in  source  properties,  attenuation,  and  site 
amplification  effects.  For  some  regions,  such  as  the  eastern  United 
States,  the  current  understanding  of  tectonic  provinces  prohibits 
localization  of  the  earthquake  source  to  particular  structures.  Hence 
the  parameters  source  size  and  distance  cannot  be  assessed  in  advance. 
Nonetheless,  the  useful  concepts  of  a  mean  and  84th  percentile  response 
spectra  have  been  developed  by  statistically  evaluating  the  accumulative 
risk  at  a  site  from  local  and  regional  sources.  Such  a  risk  analysis 
integrates  both  the  effects  of  sources  located  at  a  range  of  distances 
and  the  frequency  of  event  occurrence  for  each  magnitude  interval,  up 
to  some  maximum  source  size. 

For  the  current  project  the  concepts  of  a  mean  and  84th  percentile 
provide  a  framework  for  evaluating  the  effects  of  geologic  structure  on 
strong  ground  motion.  The  simulations  presented  below  demonstrate 
that  geologic  structure  introduces  considerable  dispersion  about  some 
mean  prediction.  For  a  particular  site  the  amplitude  bias  associated 
with  geologic  structure  may  be  either  high  or  low,  depending  upon  the 
source  location  and  the  configuration  of  the  near  site  structures.  This 
observation  suggests  that  the  problem  of  predicting  seismic  response  in 
a  basin  should  be  decomposed  into  two  components:  (1)  Assessment 
of  dispersion  associated  with  irreducible  uncertainties  (eg:  absolute 
source  location,  limited  resolution  of  geologic  structure,  uncertainties  in 
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material  properties);  and,  (2)  Systematic  site-bias,  either  high  or  low. 
This  report  focuses  on  the  second  problem.  Although  the  numerical 
results  presented  below  can  provide  valuable  input  for  the  first 
problem,  analysis  of  dispersion  must  also  include  additional  attention  to 
the  effects  of  intrinsic  attenuation  and  media  heterogeniety ,  topics  that 
are  outside  of  the  scope  of  this  project. 

Before  progressing  to  3-dimensional  modeling  of  wave  propagation  in  a 
complex  geologic  basin,  it  is  quite  useful  to  examine  the  effects  of  a 
simple  2-dimensional  model.  Figure  3.1  shows  a  single  layer  over  a 
half-space  basin  model.  The  P-wave  velocity  in  the  layer  is  2  km/sec 
and  the  basement  velocity  is  4  km/sec.  The  incident  plane  wave  energy 
is  from  the  right  side  of  the  basin,  propagating  parallel  to  the  indicated 
arrow.  Wave  propagation  through  this  model  has  been  computed  using 
a  Kirchoff  algorithm  developed  at  Sierra  Geophysics  (Apsel,  et  al, 
1983).  The  calculations  include  two  internal  multiple  reflections  within 
the  structure.  The  absolute  levels  of  motion  and  variability  of  ground 
shaking  caused  by  the  basin  structure  are  evidenced  in  the  displayed 
synthetic  seismograms  superimposed  at  seven  receiver  locations  across 
the  basin.  Note  the  long  duration  and  large  amplitudes  of  the  simulated 
records  along  the  left  half  of  the  basin  in  the  direction  of  the  incoming 
energy.  Amplification  by  the  basin  structure  caused  by  focusing  is 
responsible  for  the  amplitude  variations  of  about  a  factor  of  three. 
This  figure  clearly  illustrates  two  important  points.  Although  the 
structure  is  symetric,  a  site  on  the  edge  of  the  basin  can  be  biased 
either  high  (left  side  of  Figure  3.1)  or  low  (right  side),  depending 
upon  the  incidence  angle  and  direction  of  the  incoming  seismic  energy. 
If  the  source  location  is  unknown  in  advance,  the  most  conservative 
prediction  for  either  side  of  the  basin  should  be  based  upon  the  high 
amplitude  results  on  the  left  side  of  the  figure.  The  second  important 
effect  associated  with  the  basin  is  the  significant  increase  in  duration  of 
strong  motions  relative  to  the  flat  layered  or  rock  site  condition.  As 
the  seismic  wave  enters  into  the  basin  and  is  diffracted  (i.e.,  ray 
parameters  are  altered),  part  of  the  energy  is  effectively  trapped 
through  post-critical  reflections  at  second  and  later  bounce  points. 
This  results  in  both  a  longer  duration  of  shaking  and  higher  amplitude 
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secondary  arrivals.  The  effects  of  both  amplification  through  focusing 
and  defocusing  and  extended  duration  of  strong  shaking  through  the 
trapping  of  post-critical  reflections  are  the  two  major  effects  that 
influence  ground  motions.  The  following  paragraphs  present  the  results 
of  simulations  for  shallow  seismic  sources  located  around  a  complex 

basin  structure.  Although  the  resulting  simulated  ground  motions  show 
significant  variability  and  complexity,  the  patterns  of  motions  are  most 
easily  understood  by  referring  to  the  simple  phenomena  illustrated  in 
Figure  3.1. 

3.2  YUCCA  FLATS  -  NTS 

The  simulations  contained  in  this  section  have  been  performed  using  the 
geologic  model  developed  by  Herrin  et  al.  (1981)  for  Yucca  Flats,  NTS. 

The  principal  feature  modeled  in  this  study  has  been  the  contact 

between  the  Paleozoic  basement  rock  and  the  infilling  sediments.  Figure 
3.2  shows  the  depth  to  the  Paleozoic  contact.  This  horizon  represents 
the  most  important  impedance  contrast  affecting  amplitudes  and 
durations.  Although  this  model  will  capture  the  more  important  aspects 
of  the  simulation,  future  efforts  could  add  detail  by  incorporating 

additional  geologic  layers  within  the  basin  and  inclusion  of  the  water 
table.  Table  3.1  lists  the  physical  properties  for  each  layer  of  the 
model  used  for  the  simulation.  All  crustal  layers  below  the  Paleozoic 
contact  are  flat. 

Using  DYNARAY  complete  seismograms  were  computed  for  a  grid  of 
points  covering  the  surface  area  of  the  NTS  basin  model.  In  each 
case,  40  second  time  histories  were  developed  with  a  sample  interval  of 
0.02  seconds.  Following  the  computation  of  each  time  history,  pseudo 
velocity  response  spectra  were  computed  at  each  point  for  periods  of 
2.0  sec,  1.0  sec,  0.5  sec,  0.25  sec,  0.125  sec  and  0.0625  sec.  The 
choice  of  pseudo  velocity  response  spectra  was  chosen  specifically  to  be 
consistent  with  standard  earthquake  engineering  practices. 

Considerable  effort  was  expended  in  determining  the  ray  types  which 
produced  significant  contributions  to  the  vertical  ground  motion  in  the 
receiver  basin.  With  this  analysis,  it  was  determined  that  the  signifi- 
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cant  ground  motion  could  be  attributed  to  six  ray  types.  These  are: 

Type  I  :  The  direct  reflection  from  the  first  crustal  model 

layer  at  9  kilometers. 

Type  II  :  Same  as  Type  I  but  with  a  single  multiple  reflection 

within  the  basin. 

Type  III  :  Same  as  Type  I  but  with  a  double  multiple  reflection 

within  the  basin. 

Type  IV  :  The  direct  reflection  from  the  second  crustal  model 

layer  at  25  kilometers. 

Type  V  :  Same  as  Type  IV  but  with  a  single  multiple  reflection 

within  the  basin. 

Type  VI  :  Same  as  Type  (V  but  with  a  double  multiple  reflection 

within  the  basin. 

Figure  3.3  contains  a  schematic  diagram  of  these  six  ray  types.  All  of 
the  rays  are  entirely  compressiona!  wave  energy.  The  shear 
conversions  do  not  materially  contribute  to  the  vertical  intensities. 

The  generation  of  synthetic  time  histories  requires  both  Green's 
functions  to  model  wave  propagation  through  the  structure  and  a  source 
time  function.  Computations  with  DYNARAY  provide  both  the  requisite 
Green's  functions  and  the  convolution  of  a  user  specified  time  function 
with  each  transfer  function  to  produce  the  final  seismogram.  For  this 
study  we  have  used  an  isotropic  source  function  derived  from  numerical 
simulations  of  a  shallow  cratering  nuclear  explosion  (Trulio,  personal 
communication). 

Figures  3.4,  3.5,  and  3.6  show  the  short  period  (0.0625  sec)  response 
of  the  basin  for  three  source  locations.  The  response  values  presented 
on  these  figures  includes  geometric  attenuation  with  increasing  distance, 
changing  reflection  coefficients  of  the  crustal  reflectors  with  increasing 
distance,  and  the  site  effects  induced  by  the  basin  structure.  Each 
phenomenon  introduces  considerable  variations  in  amplitudes.  Note  that 
the  overall  response  in  Figures  3.4,  3.5,  and  3.6  generally  decay  with 
increasing  distance,  associated  with  geometric  spreading.  Superimposed 
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TABLE  3.1  PHYSICAL  PROPERTIES  OF  YUCCA  FLATS  BASIN  MODEL. 


Layer 

VP 

VS 

RHO 

ZMIN 

1 

2.400 

1.386 

2.200 

-1.172 

2 

5.000 

2.887 

2.600 

-9.000 

3 

6.300 

3.637 

3.300 

-25.000 

HS 

7.600 

4.388 

3.300 
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upon  the  decay  of  the  response  levels  with  increasing  distance  is  a  10 
km  wide  band  of  high  response  centered  at  a  distance  of  about  25  km. 
This  increase  in  ground  motion  is  associated  with  large  amplitude 
critical  reflections  from  the  first  shallow  crustal  layer  below  the 
Paleozoic  contact.  Finally,  the  structure  of  the  basin  superimposes 
additional  complexity  to  the  simulated  response. 

Figures  3.4,  3.5,  and  3.6  show  that  in  order  to  address  questions  of 
the  effects  of  shallow  basin  structure  upon  strong  ground  motion  the 
effects  of  other  major  wave  propagation  phenomena  must  be  removed 
from  the  simulated  response.  A  standard  method  for  normalizing 
response  spectra  when  investigating  the  effects  of  various  parameters  is 
to  reduce  the  response  spectra  to  a  common  distance.  The  reduction 
factor  is  the  ratio  of  the  mean  expected  values  of  the  high  frequency 
limit  of  the  response  spectra  at  the  normalizing  distance  and  the 
recording  distance.  This  procedure  references  the  response  relative  to 
the  expected  value  of  peak  ground  acceleration  (PGA).  For  application 
to  the  simulations  carried  out  in  this  study,  it  is  first  necessary  to 
develop  a  mean  expected  curve  for  PGA  for  the  NTS  -  Yucca  Flats 
basin  model.  This  has  been  accomplished  by  plotting  the  simulated  high 
frequency  response  versus  distance  for  each  source.  Next,  a  smooth 
curve  was  fit  through  the  mean  of  the  simulated  data.  This  curve 
defines  the  expected  attenuation  of  PGA  with  distance  for  the  Yucca 
Flats  model.  Finally,  the  response  spectra  computed  at  each  grid  point 
covering  the  basin,  for  each  source  location,  was  normalized  by  the 
mean  expected  PGA  to  a  common  distance.  Figures  3.7  through  3.42 
present  the  normalized  response  spectra  at  six  periods  for  six  source 
locations. 

Careful  examination  of  Figures  3.7  through  3.42  show  features  similar  to 
those  discussed  above  and  illustrated  with  Figure  3.1.  Depending  upon 
the  source  location,  many  regions  of  the  model  experience  both 
increases  and  decreases  in  the  simulated  response  spectra.  This 
results  from  the  effects  of  focusing  and  defocusing  and  will,  in  general, 
be  very  dependent  upon  the  geometry  of  the  source  and  basement 
topography.  For  instances  when  the  source  location  cannot  be  well 


SGI -R-83-096 


defined  in  advance,  one  value  of  simulations  such  as  those  presented  in 
this  report  is  to  define  the  dispersion  of  the  response  spectra  about  a 
mean  estimate.  Although  it  is  outside  of  the  scope  of  this  project, 
future  analysis  of  a  basin  response  should  include  both  a  detailed 
compilation  of  the  dispersion  and  a  comparison  of  that  dispersion  with 
observations  from  relevant  strong  ground  motion  data  sets. 

Further  examination  of  Figures  3.7  through  3.42  reveals  that  some  areas 
of  the  Yucca  Flats  basin  are  systematically  biased  either  high  or  low 
relative  to  the  mean  response.  High  amplitude  sites  are  generally 
associated  with  basement  structures  that  form  natural  "amphitheaters" 
that  focus  energy  from  most  azimuths  towards  a  small  region.  Examples 
of  this  effect  are  clearly  seen  for  the  small,  partially  closed  basin  at 
the  extreme  north  end  and  at  the  southeast  end  of  Yucca  Flats,  Figure 
3.2.  Systematic  defocusing  and  lower  average  amplitudes  appear  to  be 
associated  with  the  north-south  trending  basement  ridge  along  the 
western  side  of  the  basin,  Figure  3.2.  In  interpreting  the  simulations 
along  the  basement  ridge,  one  should  be  careful  not  to  confuse  lower 
amplitudes  from  defocusing  with  the  effective  loss  of  impedance 
amplification  due  to  basement  outcrop  over  limited  areas  along  the  crest 
of  the  ridge,  Figure  3.2. 

In  conclusion,  application  of  DYNARAY  to  the  Yucca  Flats  basin  has 
demonstrated  the  programs  capability  to  model  complex  geologic 
structures.  For  source-site  geometries  that  can  be  defined  in  advance 
and  for  well  defined  geologic  structures  this  program  can  be  used  to 
assess  site  bias.  For  the  related  problem  when  source  locations  are 
unknown  in  advance,  the  code  is  equally  useful  for  developing  the  data 
to  statistically  evaluate  both  the  dispersion  and  the  probability  that  a 
specific  site  response  may  be  biased  high  or  low. 

This  work  has  led  to  the  identification  of  several  desirable  features  of 
DYNARAY  that  would  extend  the  usefulness  of  the  code.  These  topics 
are  addressed  in  the  recommendations  section  of  this  report. 
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4.0  RECOMMENDATIONS 

This  project  has  resulted  in  the  very  successful  development  of  both  a 
theoretical  basis  for  dynamic  raytracing  and  a  functional  computer 
program,  DYNARAY.  During  the  application  of  DYNARAY  to  the 
modeling  of  strong  ground  motion  at  Yucca  Flats  (Section  3  of  this 
report)  several  recommendations  for  future  enhancements  were 
developed  that  would  improve  the  usefulness  of  DYNARAY  beyond  that 
envisioned  in  the  original  statement  of  work.  These  enhancements  are 
outlined  below. 

4.1  INCORPORATION  OF  A  VELOCITY  GRADIENT  IN  THE  HALF-SPACE 
For  very  shallow  seismic  sources,  such  as  those  modeled  in  this  study, 
the  principal  energy  that  reaches  a  site  within  a  large  basin  comes  from 
raypaths  that  are  refracted  or  turned  around  in  the  shallow  crust. 
Using  DYNARAY,  these  shallow  raypaths  can  be  modeled  with  a  series 
of  layers  that  simulate  a  velocity  gradient.  However,  if  the  layers  are 
relatively  coarse,  then  the  expected  smooth  amplitude  decay  with 
distance  will  be  irregular  owing  to  the  location  of  critical  reflection 
angles.  A  valuable  addition  to  DYNARAY  would  be  the  inclusion  of  a 
linear  velocity  gradient  in  the  model  half-space.  This  would  provide  a 
capability  for  inexpensively  modeling  turning  rays  and  would  improve 
the  accuracy  of  the  simulated  strong  ground  motion.  Furthermore,  this 
capability  would  provide  a  means  to  completely  separate  the  effects  of 
basin  structural  amplification  from  the  mid-crustal  reflections. 

4.2  PHYSICAL  ATTENUATION 

Elastic  body  waves  decay  as  R  \  However,  it  is  well  recognized  that 
strong  ground  motion  parameters  such  as  peak  ground  acceleration 
decay  as  R  to  R  .  This  discrepency  translates  into  substantial 
differences  in  predicted  motions  at  distances  beyond  a  few  kilometers. 
The  difference  in  decay  rates  between  the  elastic  case  and  the  real 
world  is  well  modeled  by  the  addition  of  physical  attenuation  (see  for 
example  Hadley  and  others,  1982).  As  DYNARAY  computes  travel  time 
within  each  layer  for  each  ray,  a  final  attenuation  operation  (t*)  could 
be  constructed  from  the  average  attenuation  along  each  raypath.  The 
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addition  of  this  operation  would  greatly  improve  the  accuracy  of  the 
simulated  strong  ground  motions. 

4.3  COMPILATION  OF  DISPERSION  STATISTICS 

One  of  the  valuable  parameters  for  seismic  design  that  can  be  assessed 
with  DYNARAY  is  a  measure  of  uncertainty  or  dispersion  (see  Section 
3.1).  Although  DYNARAY  outputs  disk  files  that  can  be  sorted  and 
the  required  mean  and  84th  percentile  information  extracted,  an 
additional  module  in  DYNARAY  could  efficiently  provide  this  data  to  an 
analyst. 

4.4  COMPARISON  WITH  FIELD  OBSERVATIONS 

Before  DYNARAY  is  used  extensively  for  assessing  the  effects  of 
geologic  structure  on  strong  ground  motion,  a  limited  field  test  of  the 
program  would  be  valuable.  This  would  provide  both  great  insight  into 
the  degree  of  geologic  detail  required  to  adequately  model  the  effects  of 
structure  and  provide  a  strong  validation  of  the  accuracy  of  the  code. 
The  field  experiment  envisioned  would  involve  the  deployment  of  field 
recording  digital  siesmographs  across  a  known  geologic  basin.  Next 
small  shots  would  be  fired  at  points  both  within  and  around  the  basin. 
The  resulting  seismic  records  would  provide  an  ideal  dataset  for 
comparison  with  the  DYNARAY  simulations. 
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Figure  3.1  Waveforms  and  relative  amplitudes  for  a  wave  incident  at  a 
shallow  angle  from  the  right  (arrow).  The  synthetic  time 
histories  have  been  calculated  with  the  Kirchhoff  technique 
and  includes  two  multiples  within  the  basin.  Note  the 
factor  of  3  increase  in  amplitude  caused  by  structure. 
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Figure  3.3  Raypaths  used  in  the  Yucca  Flats  Study 
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Figure  3.4  Variations  in  short  period  response  across  Yucca 
Flats  for  source  location  shown  by  symbol  "S". 
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Figure  3.5 


Variations  in  short  period  response  across  Yucca 
Flats  for  source  location  shown  by  symbol  "S". 
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Figure  3.6  Variations  in  short  period  response  across  Yucca 
Flats  for  source  location  shown  by  symbol  "S". 
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Figure  3.7  Variations  in  pseudo  velocity  response  across 

Yucca  Flats  for  the  period  shown.  The  source 
location  is  indicated  by  the  symbol  "S".  As 
discussed  in  the  text,  the  response  has  been 
normalized  by  the  mean  expected  PGA  for  the 
Yucca  Flats  model. 
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Figure  3.8  See  Figure  3.7 
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Figure  3.10 


See  Figure  3.7 
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YUCCA  FLATS  BASIN  MODEL 
PERIOD  =  0.  1250  SECONDS 
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Figure  3.12  See  Figure  3.7 
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Figure  3.13 


See  Figure  3.7 
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See  Figure  3.7 
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YUCCA  FLATS  BASIN  MODEL 
PERIOD  =  0.  2500  SECONDS 
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Figure  3.16  See  Figure  3.7 
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YUCCA  FLATS  BASIN  MODEL 
PERIOD  =  0. 1250  SECONDS 


EAST  -  WEST  DISTANCE  IN  KILOMETERS 
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Figure  3.17 


See  Figure  3.7 
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YUCCA  FLATS  BASIN  MODEL 
PERIOD  =  0.  0625  SECONDS 

EAST  -  WEST  DISTANCE  IN  KILOMETERS 
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Figure  3.18 


See  Figure  3.7 


zooct-i  -  inooi-r  cu-*ini-<tzuuj  <-<z  acw-iDziuh-uctn 


YUCCA  FLATS  BASIN  MODEL 
PERIOD  =  2.  0000  SECONDS 


53 


EAST  -  WEST  DISTANCE  IN  KILOMETERS 


10. 


20. 


30. 


40. 


0.  10.  20  30 

1.  O 

#### 
#### 
#### 
#### 

O.  9 

7X7.7. 
7.7.7.'/. 
7X7.7. 
7.7X7. 

O.  8 

O.  7 

XXXX 
XXXX 
XXXX 
XXXX 

0  6 

xxxx 

XXXX 
XXXX 

xxxx 

O.  0 

++++ 
++++ 
++++ 
++++ 

O.  4 

\\\\ 
\\\\ 
\\\\ 
\\\\ 

O.  3 

//// 
//// 
//// 
//// 

O.  2 


VALUE  =  SCAL*  O. 327E+01  +  O. 970E-02 


/// 

:  :  / 


:  :  .  :  : 

:  .  .  /  \\\/ 

: : \/  \\/ 

:  /:  :  \/.  ////\ 

//: : /////\/ 

/  .  .  /////\\  \ 

/  .  .  ///++\/ 

/.  :  :  ///\\\/ 

\/ :  :  :  ////\/ 

//: : : //\\\\x  / 

\/. . . //\\/\\// 

#/. . //////\/\\/ 
///+/\\\///\\/ 
///: //\\\\///// 
////: ////\///\/ 
/////////\///\/ 
////: ///\\//\\ 
///  : \\/\\//\: \\ 
///  :  W//WW:  .  . 
\\//. ///////\: 

//: ///////// 

: : ///////// 

: : //////// 

: : //////// 

: : ////////: 

. . \: : /////\ 

.  ///\\/ 

.  :  .  .  .  //////.  . 

.  :  :  ////.  . 

.  :  :  ///:  .  . 


-++ 


O.  1 


O.  0 


Figure  3.19 


See  Figure  3.7 
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Figure  3.20 


See  Figure  3.7 
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Figure  3.21 


See  Figure  3.7 


SGI-R-83-096 


cnam-inirsarn?:  2m  rmzJHKi«D  ncoui  -  naoz 


YUCCA  FLATS  BASIN  MODEL 
PERIOD  =  0.2500  SECONDS 


56 


EAST  -  WEST  DISTANCE  IN  KILOMETERS 


lO. 


20. 


30. 


40. 


0. 

10 

20. 

30. 

+ _ 

1.  0 

• 

■ 

1 

#### 

• 

i 

1 

4444 

« 

• 

1 

4444 

« 

• 

t 

4444 

• 

< 

//  .  :  :  V 

« 

1 

0.  9 

:  /  \\\/ 

t 

t 

Y.7.7.7. 

« 

1 

.  \/  \\/ 

» 

1 

V.7.7.7. 

:  :  /// 

:  \/.  ////\ 

• 

1 

7.V.7.7. 

t 

• 

.  \////\/ 

• 

1 

7.V.7.7. 

.  .  :  / 

.  ///\\\/  \ 

+ 

O.  8 

1 

1 

:  :  +/ 

/////\/ 

1 

*#** 

< 

:  :  //. 

:  /////\/ 

* 

1 

**** 

• 

1 

. //: \/ 

:  :  ////\/ 

i 

1 

**** 

« 

1 

///// 

:  :  //\\\\x  \ 

i 

i 

**#* 

1 

« 

:  :  :  \/ 

.  .  //\\/+\//  S 

i 

1 

0.  7 

« 

.  #/ 

.  //////\/\\/ 

• 

1 

xxxx 

« 

/ 

:  +/\\\///\\/ 

* 

1 

xxxx 

1 

. ///: //\\\\/\\// 

* 

1 

xxxx 

1 

1 

;  : 

///:  ////\/\\\/ 

» 

1 

xxxx 

+ 

////\ : ////\///\/ 

0.  6 

« 

1 

////  ///\\//\\ 

• 

1 

xxxx 

1 

1 

:  /: 

///  :  \\/\\//\. \\ 

• 

l 

xxxx 

« 

1 

:  /. 

///  .  \\//\//+;  .  . 

i 

1 

xxxx 

1 

1 

\\//: ///////\: 

• 

1 

xxxx 

1 

1 

: : //: ///////// 

i 

1 

O.  3 

1 

1 

//:  :  :  /■//////// 

i 

1 

++++ 

• 

1 

:  :  ////////' 

» 

1 

+  +  +  + 

1 

: : //////// 

» 

1 

+  +  +  + 

• 

1 

.  .  :  :  ////////: 

i 

l 

+ 

.  .  .  .  \: : f ff  f/\ 

•f* 

O.  4 

« 

1 

.  ////// 

1 

\\\\ 

• 

< 

.  .  .  :  .  .  .  //////.  . 

• 

I 

\\\\ 

1 

1 

.  :  :  ////.  . 

i 

i 

\\\\ 

» 

« 

i 

\\\\ 

1 

« 

l 

O.  3 

1 

i 

i 

//// 

t 

1 

• 

1 

//// 

■ 

• 

1 

/  f  f  r 

1 

1 

» 

1 

//// 

+ - + - + - ++  O  2 


VALUE  =  SCAL*  O.  327E+02  O.  984E-01 


O.  1 


O.  O 

Figure  3.22  See  Figure  3.7 
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Figure  3.23 


See  Figure  3.7 
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Figure  3.24 


See  Figure  3.7 


SGI  - R -83-096 


zoki-x  --  wooe-x  q  ■->  in  t- «  z  u  m  <-* z  ac<-* joziut-woEtn 


59 


YUCCA  FLATS  BASIN  MODEL 
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Figure  3.25 


See  Figure  3.7 
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Figure  3  .  ?b  See  Figure  3.7 
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Figure  3.27  See  Figure  3.7 
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Figure  3.28  See  Figure  3.7 
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See  Figure  3.7 
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Figure  3.31  See  Figure  3.7 
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Figure  3.32  See  Figure  3.7 
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Figure  3.33  See  Figure  3.7 


SGI-R-83-096 


z>->  —  a  ucacn  —  i-cnaz 


VUCCA  FLATS  BASIN  MODEL 
PERIOD  =  0.  2500  SECONDS 


EAST  -  WEST  DISTANCE  IN  KILOMETERS 


0. 

10. 

20. 

30. 

•  ! 

•tut  \ 

\  ### 

\\. 

+  +X  +  +\ 
\\+: /\\ : : + 
:  \  .  /  \  \  \ : 

.  f.  77*. 


\ \ \ \ / : . 
\++\\\//: /  . 
+x\\'//: . . 

/  .  .  -*++ \\\\\// 

+  //xx  X++N//W. 
+//\\+xx+\\\\\ 
/\///\\\\++\\\+/ 
/\//\\\\\\++xx+. 
\\\//\\\\\\\W\X  + 
\+\++/  \++\\\\\\X++ 

.  x \ x  x /  \++\\\\\+. 
++++/\\\\\\\\\\/ 
+++\\+++\\\W/ 

X  x++-*-+++\\\\\  + 
++x++++++\\\/: 
x++x  x  ++\ \\7 .  .  .  . 
++x  x+++\\\:  .  . 
+xxxx++++V 
X++X++++V  .  . 
++xX++++\ 
x  x  x  x  x+++\\\: 
xx  Xx+++\++. 
.  . .  +\++. / 


#### 

#### 

#### 

#### 

7.7.7.  V. 

77.7.7. 
77.7.7. 
T.7.7.7. 

**** 

**■»* 

xxxx 

xxxx 

xxxx 

xxxx 


++++ 

++++ 

\\\\ 

\\\\ 

\\\\ 

\\\\ 

7  /  7  / 
7  7  7  7 
7  7  7  7 
7  7  7  7 


-++  O.  2 


VALUE  =  SCAL*  O.  184E+02  +  O.  127E+00 


Figure  3.34 


See  Figure  3.7 
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Figure  3.35  See  Figure  3.7 
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See  F igure  3 . 7 
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